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I. INTRODUCTION 

The atomic vibrations that are always present in 
condensed matter lead to important corrections to the 
static-lattice model of crystalline solids. The most com- 
mon approach for evaluating such corrections is to in- 
clude vibrational effects within the harmonic phonon 
approximation,^^ which works well in many circum- 
stances. However, anharmonic effects can be impor- 
tant in, for example, systems incorporating light atoms 
which have large vibrational amplitudes, systems with 
weak bonding such as hydrogen bonding, and systems 
at high temperatures. When the vibrations are present, 
their effect on the electrons can also become important. 
Electron-phonon interactions are central, for example, to 
the description of superconductivity. 

Renormalization of the phonon degrees of freedom is an 
important idea underlying most treatments of phonon- 
phonon and electron-phonon interactions. For example: 
(i) the phonon frequencies as a function of temperature 
may be used to investigate the stability of crystals with 
unstable harmonic phonons, (ii) the change in phonon 
frequencies due to electron-phonon interactions may be 
used in the study of the temperature dependence of elec- 
tronic band gaps, and (iii) the phonon frequencies at dif- 
ferent volumes may be used within the quasiharmonic 
approximation to study thermal expansion. 

The self-consistent phonon method was devised by 
Born and HootonJ^l and extended by ChoquarcP and 
WerthamerP In the self-consistent phonon method, the 
phonon frequencies are renormalized by considering the 
atomic vibrations over a region about equilibrium which 
can include anharmonicity. For example, it is useful in 
systems with a harmonic instability, such as solid helium 
which forms under pressure, where including anharmonic 
terms leads to the correct physical picture in which the 
body-centered-cubic structure is stabilized. A historical 
overview of these developments is given by Klein and 
HortonP 



More recently, these ideas have been used in first- 
principles density functional theory (DFT) electronic 
structure calculations and have shown the occurrence 
of dynamically unstable phases at the harmonic level 
that are stabilized at finite temperature by means of 
anharmonicity.^ Souvatzis and co-workers P and Hell- 
man and co-workers^ have calculated renormalized 
phonons at finite temperature for systems with harmonic 
phonon instabilities. Antolin and co-workers^ presented 
a simpler and computationally less demanding method 
for solving this problem, using a harmonic fit to large- 
amplitude phonon displacements. 

Thermal expansion arises from anharmonicity 
and is usually addressed within the quasiharmonic 
approximation,^ i n which the phonon frequencies are 
taken to depend on the volume. The minimization 
of the free energy with respect to volume at a given 
temperature leads to the equilibrium volume at that 
temperature. This method has been very successful, 
but for systems with several lattice parameters it is 
computationally expensive to relax the structures at 
each volume. 

Lattice vibrations also play a central role in electron- 
phonon interactions. The study of electron-phonon in- 
teractions within many-body perturbation theory is well 
established.^ The effects of electron-phonon interactions 
on the temperature dependence of band gaps are sub- 
stantial, and they have been described within the har- 
monic approximationP^^ Recent theoretical advances 
have enabled such calculations to be performed from first 
principles.^ 

In this paper we describe a unified approach for study- 
ing the effects of anharmonicity at zero and finite tem- 
peratures on quantities including phase stability, the to- 
tal energy, electron-phonon coupling leading to band gap 
renormalization, and thermal expansion. 

The starting point of our method involves perform- 
ing first-principles calculations of the low-energy portion 
of the Born-Oppenheimer (BO)P^ energy surface defined 
by the harmonic phonon coordinates, and exploring it 
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out to large atomic displacements where the harmonic 
approximation is no longer accurate. This information 
is used to calculate anharmonic phonon free energies 
and the associated vibrational excited states using the 
vibrational self-consistent-field (VSCF) equations.^ We 
then use aperturbation theory constructed on the VSCF 
equations^ to obtain a second-order correction to the to- 
tal vibrational free energy. This m ethod has been used 
before in studying molecule d 22 * 24 1 and extended model 
systems,^ but, as far as we are aware, this is the first 
first-principles application to three-dimensional periodic 
solids. 

We have used the anharmonic wave functions to calcu- 
late quantum-mechanical expectation values of phonon- 
dependent quantities at zero and finite temperatures. For 
example, renormalization of the electronic eigenvalues 
due to their interaction with the phonons leads to the 
temperature dependence of the electronic band structure. 
Similarly, the vibrational stress tensor can be calculated 
and used to directly study thermal expansion. 

The main computational cost in our approach is that 
of performing the electronic structure calculations to ex- 
plore the BO energy surface. While performing these 
calculations, data for all quantities of interest can be ac- 
cumulated, and the further analysis of particular quan- 
tities can be done with little additional computational 
cost. 

We report results for the anharmonic energy, the 
temperature dependence of the thermal band gap, and 
the linear coefficient of thermal expansion of diamond, 
lithium hydride (H 7 Li) and lithium deuteride (D 7 Li). 
Our results are in good agreement with experiment and 
other theoretical approaches. As far as we are aware, 
no experimental or theoretical results for the tempera- 
ture dependence of the band gap of lithium hydride have 
been published previously. 

The remainder of this paper is arranged as follows. 
In Sec. ITT] we present the theoretical framework used to 



study anharmonicity and in Sec. Ill we extend the for- 
malism to the evaluation of general phonon expectation 
values. In Sec. |IV| we describe our computational ap- 
proach. Our results are presented in Sec. [V] and we draw 
our conclusions in Sec. VI All equations are given in 



Hartree atomic units, in which the Dirac constant, the 
electronic charge and mass, and An times the permittiv- 
ity of free space are unity (h = |e| = m c = Aire® = !)• 



II. ANHARMONIC TOTAL ENERGY 



ing to flei(R)|tf(R)) = £ei(R)|tf(R)). In the following 
we consider only the electronic ground state. 

The nuclear or vibrational Hamiltonian H v n, is 
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where the R p are the position vectors of the unit cells 
that make up the periodic supercell, a labels the different 
atoms within a unit cell, and m a is the mass of atom a. 
In this context, the electronic eigenvalue as a function of 
atomic position £7 e i(R) is called the BO surface. 



B. Harmonic approximation 

Within the harmonic approximation (HA ^3 the BO 
energy surface is approximated as a quadratic function 



of the atomic displacement coordinates u„ Q = r. 
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where r pa are the atomic positions, r pa the equilibrium 
atomic positions, and Latin indices i and j are used to 
label Cartesian directions. The HA often works very well 
because the nuclei are heavy and therefore only explore 
the neighborhood of their equilibrium positions. In the 
following we drop the constant energy term E e \(R°). 

The matrix of force constants is defined as 
Aa;ja'(Rp,Rp) = d 2 E c \ (R°)/<9 u pQ;i d u p : a >- 3 , , and the 
dynamical matrix at the reciprocal space point k is 



Dia\ja' (k) — 
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where N p is the number of unit cells within the periodic 
supercell. 

The vibrational Hamiltonian of Eq. ([2]) can be rewrit- 
ten in terms of normal or phonon coordinates q n k, which 
are related to displacement coordinates by 



A. BO approximation 

Within the BO approximation^!! the electronic and 
vibrational degrees of freedom are treated separately. 
The atomic positions parameterize a family of electronic 
Hamiltonians H C \(H) labeled by the collective nuclear po- 
sition variable R. Each of these Hamiltonians leads to 
the electronic energy E c \ and wave function |\&) accord- 
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Here, Wk n -i a are the eigenvectors of the dynamical matrix 
and n is the phonon branch index. In terms of phonon 
coordinates, which can always be chosen to be realj^ the 
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vibrational Hamiltonian is 
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where are the eigenvalues of the dynamical ma- 
trix. This Hamiltonian consists of a sum of terms for 
non-interacting simple harmonic oscillators of frequen- 
cies W„k- 



C. Anharmonic approximation 

1. Principal axes approximation to the BO energy surface 

To study anharmonic properties it is necessary to find 
an approximation to the BO energy surface that goes 
beyond the HA. The HA is often very accurate, and in 
such cases anharmonicity is expected to be a perturba- 
tion about the harmonic case. We write the BO energy 
surface as^ 

E cl (Q) = E cl (0) + J2 K*k(<7„k) 

n , k 

+ Kik;n'k'(<?nk)QWk') H , (7) 

n,k n',k' 

where Q is a collective phonon vector with elements q n k, 
the primed sum indicates that the term (n, k) = (n' , k') is 
excluded and the factor 1/2 accounts for double counting. 
The independent phonon term takes the form 



V„k(<?nk) = E c \(0, . . . , q„ k , • ■ • , 0) - B e i(0) 



(8) 



but note that we do not assume the HA. The two-body 
term which introduces coupling between the modes takes 
the form 

Kik;n'k'(<Zrak, <7n'k' ) =E e i(0, . . . , q n k 7 . . . , q n 'k' , ■ ■ ■ , 0) 

-14k(9nk)-K'k'(9n'k')-^el(0). (9) 

It is possible to continue this series for more general N- 
body terms, but because we start from the HA in which 
phonons are exactly noninteracting, it is expected that 
the higher-order terms will decrease in magnitude as N 
increases. [We note that Eq. (14) of Ref. [M] contains a 
typographical error connected with the coupling of vibra- 
tional modes.] 

We show examples of the BO energy surface from inde- 
pendent phonon terms in Fig. [I] for diamond and lithium 
hydride. Both plots show the BO energy surface as 
a function of phonon amplitude for an optical phonon 
at the center of the Brillouin zone. The potential for 
diamond is asymmetric and the dominant anharmonic 
term is cubic. Diamond has a face-centered cubic (fee) 
crystal structure with a two-atom basis at (0, 0, 0) and 
(a/4, a/4, a/4), and a conventional cubic cell of lattice pa- 
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FIG. 1. (Color online) Harmonic and anharmonic BO energy 
surfaces per unit cell for an optical F-point phonon in (a) 
diamond and (b) H 7 Li. The leading anharmonic term for di- 
amond is cubic, and that for H 7 Li is quartic. For a harmonic 
ground-state wave function (Gaussian), the one-sigma ampli- 
tude is 9.1 a.u. for diamond and 13.7 a.u. for H 7 Li. In real 
space, the one-sigma amplitudes correspond to 0.043 a.u. for 
diamond and 0.299 a.u. for H 7 Li. 



rameter a. The optical phonon represents the motion of 
the two atoms of the basis in opposite directions, leading 
to a steep potential when they approach, and a shallower 
potential when they move apart. The lithium hydride 
crystal structure is also fee with a two-atom basis, but 
the second atom is located at (a/2, a/2, a/2), leading to 
a rather isotropic environment and a symmetric anhar- 
monic potential. In this case the anharmonicity arises 
because the potential is steeper than the harmonic one, 
and the leading anharmonic term is quartic. 

The crystal structure determines the symmetry opera- 
tions relating different k-points within the first Brillouin 
zone. This implies that the only points that need to be 
treated explicitly are those within the irreducible wedge 
of the Brillouin zone, which reduces the computational 
cost significantly. For example, a 5 x 5 x 5 supercell of 
the diamond structure has 125 k-points, of which only 22 
need be treated explicitly. 



2. Vibrational self- consistent field 

Within the principal axes approximation to the BO 
energy surface the vibrational Hamiltonian is 



#vib = 
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(10) 



and the Schrodinger equation is solved by an anharmonic 
vibrational wave function |$(Q)) and corresponding an- 
harmonic energy E vib , where H vih \$(Q)) = E vih \®(Q)). 
In Eq. (10), /? is the inverse temperature. In principle it 



would be straightforward to calculate the BO energy sur- 
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face with electrons distributed, for example, according to 
the Fermi-Dirac distribution at some finite temperature 
without substantial additional cost. This temperature 
dependence can be important for metallic systems, but 
we will drop it in the rest of this paper, because here we 
only study systems with a band gap. 

To solve this equation, the energy is minimized with 
respect to a set of single-phonon states {|</> n k(<7nk))} 
that form a Hartree product for the trial wave function 
l^(Q)) = Iln k l^nk^nk))- This leads to the vibrational 
self-consistent field equations^! 
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where A„k is a vibrational energy eigenvalue and 
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The primed products indicate that the term (n, k) is ex- 
cluded. The final energy is 
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The VSCF equations determine a set of excited eigen- 
states for each degree of freedom. These can be used 
to construct approximate anharmonic excited state wave 
functions as 



i* s (Q)>-ni^ k (^)), 



(14) 
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where the superindices 5* n k label the excited states, and 
S is a vector with elements S n k. The corresponding en- 
ergies are labeled by Es- 

A perturbation theory can be constructed on the VSCF 
equations^, and the second order correction to the en- 
ergy of state S is 
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The anharmonic free energy F can be calculated at any 
inverse temperature f3 from the eigenvalues that solve the 
VSCF equations as 



F 
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InZ, 
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where Z — ^ s e @ Es is the partition function. 



III. PHONON EXPECTATION VALUES 

A. General formulation 

The wave function |$(Q)) that solves the VSCF equa- 
tions contains, in principle, all of the physical informa- 
tion about the vibrational system within the principal 
axes approximation. For an operator O(Q) that depends 
on the phonon coordinates, it is then possible to evaluate 
the expectation value with respect to the phonon wave 
function as 



(0(Q)>* = ($(Q)|6(Q)|$(Q)>. 



(17) 



The expansion of the operator O(Q) in terms of the 
phonon amplitudes can be written in analogy to the ex- 
pression for the BO energy surface, as 



O(Q) = 6(0)+£6 nk (g nk ) 



^ ^ 0„k;n'k'( ( J , nk, 9n'k') 
,k n'.k' 



(18) 



This expression can be constructed from data accumu- 
lated during the mapping of the BO energy surface along 
the principal axes. 

Within the finite temperature formalism, the expecta- 
tion value of the operator O(Q) at inverse temperature 
P is 



<0(Q))<M = \ $> s (Q)|0(Q)|<i> s (Q)) e 



(19) 



Examples of phonon-dependent operators are elec- 
tronic eigenvalues, the stress tensor, and the atomic po- 
sitions. In this paper we will describe electronic eigen- 
values and their relation to electron-phonon interactions 
and the stress tensor and its relation to thermal expan- 



B. Electron-phonon interactions 

The effects of electron-phonon interactions on the elec- 
tronic band structure of solids can be calculated by renor- 
malizing the phonons due to interactions with the elec- 
tronic states^SED or renormalizing the electrons due to 
their interactions with the phononsPSl These two ap- 
proaches can be shown to be equivalent by Brooks' 
theoremPS 

We accumulate the single-particle electronic band 
structure while mapping the BO energy surface within 
the principal axes approximation. Using the expansion 
above for the phonon amplitude dependence, it is then 
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possible to calculate the renormalized electronic band 
structure due to the presence of the phonons. 



Stress tensor 



The differential Gibbs free ener gy d G of a system with 



an externally applied stress erf 



dG = dF cl + dF vih Vi? de i 



(20) 



where is the strain tensor and fi is the volume. The 
internal vibrational stress tensor is defined as 



1 dF vi 



OV; = - 



vib 



(21) 



The vibrational energy is an approximately linear func- 
tion of volume for the configurations of interest in ther- 
mal expansion. This leads to a constant vibrational stress 
tensor. It is then possible to define an effective stress 
cr? ff = crff* + ajf 3 such that Eq. (po| can be rewritten as 



(22) 



The minimum of the Gibbs free energy with respect to 
variations in strain corresponds to the equilibrium config- 
uration of the system. Therefore, minimizing the Gibbs 
free energy in Eq. ( 22 ), where the vibrational dependence 
is implicit in the effective stress erf?, leads to the equi- 
librium configuration for the full vibrational system. If 
the vibrational stress varies significantly in the region of 
interest, an iterative approach can be adopted. Equation 



( 22 1 can be used to determine a new equilibrium config- 



uration, which can be used as a starting point for the 
next calculation. This iterative process can be repeated 
until convergence is reached. In this work, only one or 
two iterations were required for convergence. 

The stress tensor arising from the electronic Hamilto- 
nian can be calculated for each point sampled on the BO 
energy surface, and its expectation value with respect to 
the phonon wave function evaluated using Eq. ( 19 ). This 



stress, cry ' , is the contribution of the potential energy 
to the vibrational stress. In general, the internal stress 
tensor has contributions from both the kinetic and po- 
tential parts of the Hamiltonian. The contribution from 
the kinetic partp^ 
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is therefore also required. The symbol u indicates the 
time derivative of a displacement coordinate. For an 
isotropic system where the internal vibrational stress can 
be described by a pressure P, the contribution from the 



kinetic part becomes 3P = — tr er 



vib.T 



2Ej ih /n, where 



tr denotes the trace of the tensor and Ej ih is the vibra- 



vib 



ib,T 



tional kinetic energy. 

The total vibrational stress a 
then be used in tandem with Eq. (|22|) to determine the 



vib,V 
r . . 

y 



equilibrium configuration in the presence of vibrations 
and at finite temperature. 



IV. COMPUTATIONAL DETAILS 

A. Electronic calculations 

We have solved the electro nic S chrodinger equa- 
tion within plane-wave using ultrasoft 
pseudopot ent ialiP^ as implemented in the CASTEP 
codePS All energy differences were converged to within 
10~ 4 eV per unit cell and all stresses were converged to 
within 1CP 2 GPa. This level of convergence required 
an energy cutoff of 1000 eV and a reciprocal-space 
Monkhorst-PaclPSl grid f spacing 2tt x 0.04 A -1 We 
have used the local density approximation (LDAp£122l 
to the exchange-correlation functional for the diamond 
calculations, and the Perdew-Burke-Ernzerhof (PBE) 38 
generalized gradient approximation density functional 
for the lithium hydride and dcutcridc calculations. 



B. Vibrational calculations 

We have mapped the BO energy surface along the 
phonon modes with a maximum amplitude given by a 
multiple (between 3 and 5) of the harmonic expectation 
value of the mode amplitude 
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(24) 



We have fitted the independent phonon term with a poly- 
nomial of order 6, and tests with polynomials of or- 
der up to 8 have confirmed the convergence of the re- 
sults. For the two-body terms we have fitted the two- 
parameter functional form V n wk> = c ik ; „'k' 'JVikgn'k' + 

c ikn'k' ^nk^n'k'- The nrst term is the lowest order cou- 
pling term in a Taylor expansion of the potential. The 
second term is not the next order term, but the low- 
est order term that gives a non-zero contribution when 
overlapped with the harmonic ground state, and, for this 
reason, it is expected to be larger than other lower order 
cross-terms. Numerical tests have confirmed this. 

We use simple harmonic oscillator eigenstates as a ba- 
sis set for expanding the anharmonic wave function in 
the VSCF equations. The simple harmonic oscillator fre- 
quencies used are those of a quadratic fit to the anhar- 
monic independent-phonon BO energy surface. We have 
used basis sets including 100 states for each mode, check- 
ing the convergence of the results by using larger sets. 
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FIG. 2. (Color online) Anharmonic energy correction as a 
function of the number of normal modes (supercell size) for 
diamond, H 7 Li and D 7 Li. The anharmonic correction is larger 
for the lighter elements. 



V. RESULTS 
A. Anharmonic energy 

The mapping of the BO energy surface within the prin- 
cipal axes approximation requires the displacement of the 
atomic nuclei inside a periodic supercell. The supercell 
size determines a set of commensurate reciprocal-space 
points associated with the harmonic phonon modes. Con- 
vergence with respect to supercell size must be tested. 

The sampling of the Brillouin zone within the HA is 
not constrained by the supercell size in the same man- 
ner. This is because within the HA, all that is needed is 
the dynamical matrix at the k-points of interest, and this 
matrix can be constructed at general k-points irrespec- 
tive of the supercell size. This is accomplished by first 
constructing the real-space matrix of force constants and 
then Fourier transforming to the dynamical matrix at a 
general k-point. This approach relies on the fact that 
the force constants tend rapidly to zero as the separation 
between atoms increases. The matrix of force constants 
can be constructed using a finite displacement method 39 
or density functional perturbation theorj^UiU to evalu- 
ate the dynamical matrices on a coarse k-point grid, and 
then calculating the matrix of force constants by means 
of an inverse Fourier transform. 

The strategy we follow is to calculate both the har- 
monic -Fhar an d anharmonic -Fanh free energies from the 
k-points commensurate with the supercell. The anhar- 
monic correction is then evaluated as AF an h = F an h — 
.Fhar, and converged with respect to supercell size. This 
correction is then added to an independently converged 
harmonic free energy to give the total anharmonic energy. 

The supercell-size dependence of AF an h at zero tem- 
perature for diamond, lithium hydride, and lithium deu- 
teride is shown in Fig. [2] The anharmonic energy is eval- 



uated using only the independent phonon terms in Eq. 

The energy difference between the last two succes- 
sive points in each curve is smaller than 0.2 meV per unit 
cell. Table |T] gives the numerical values of the harmonic 
phonon energy and the independent phonon anharmonic 
correction per unit cell. 

The anharmonicity of diamond is small. For the small- 
est supercell size corresponding to the T-point only, the 
anharmonic correction is negative, as expected for cubic 
anharmonicitjEl ( see pjg nj Including more k-points 
reverses the sign of the anharmonic correction, but it re- 
mains small. 

The different isotopic masses of H 7 Li and D 7 Li lead to 
different vibrational properties. The anharmonic correc- 
tion is larger for H 7 Li than for D 7 Li because the lighter 
elements explore wider regions of the BO energy surface 
and encounter more anharmonicity. We have also per- 
formed calculations with the lithium isotope 6 Li, but the 
results are not much affected by this isotopic substitu- 
tion because it represents a mass reduction of only about 
15%, compared to the doubling of mass from substituting 
hydrogen by deuterium. 

We have studied the effects of including the phonon 
coupling term in Eq. ^ for a 2 x 2 x 2 supercell (48 
modes) of H 7 Li. The anharmonic coupling correction 
has the opposite effect to the independent phonon cor- 
rection, and the final anharmonic correction with cou- 
pling between modes is AF an h = +0.79 meV, compared 
to AFanh = +3.21 meV with independent phonons only 
(see Table M. U sing the second order perturbation the- 
ory of Eq. (|15| does not change the final result appre- 
ciably, indicating that the energies have converged with 
respect to the mean-field theory. This means that, like 
diamond, H 7 Li has a small anharmonic energy correc- 
tion, but, unlike diamond, the reason for this is the can- 
cellation of the contributions from independent phonons 
and pairwise phonon coupling. This negative contribu- 
tion from the coupling terms can be better appreciated 
by referring to Fig. [3] We show the coupling BO en- 
ergy surface Vi-2 (91,92) for two optical phonons (91,92) 
at the T-point compared to the approximate surface con- 
structed from the independent phonon terms only. The 
accurate potential is shallower than the approximate one, 
leading to the smaller anharmonic correction when in- 
cluding phonon coupling. The corresponding harmonic 
two-dimensional subspace is very similar to the approx- 



TABLE I. Converged anharmonic free energy correction per 
unit cell at zero temperature using independent phonons for 
the mapping of the BO energy surface along the principal 
axes. The converged harmonic energy is also displayed. 





Num. modes 


F har (meV) 


AF anh (meV) 


Diamond 


750 


367.59 


0.23 


H 7 Li 


384 


221.49 


3.21 


D 7 Li 


384 


174.75 


1.88 
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FIG. 3. (Color online) Anharmonic two-dimensional BO en- 
ergy surface per unit cell for a pair (qi, (72) of optical modes at 
the T-point for H 7 Li. We plot (a) the exact energy surface and 
(b) the approximate energy surface for independent phonons 
only. For a harmonic ground-state wave function (Gaussian), 
the one-sigma amplitude is 13.7 a.u. 



imate (bottom) potential in Fig. |3j but with concentric 
circular equipotcntial lines because the harmonic phonon 
frequencies are degenerate in the case shown. 

Nolan and co-workers^ have used very accurate elec- 
tronic structure methods to calculate the cohesive energy 
of lithium hydride, including the zero point (ZP) energy 
at the harmonic level. Their accuracy is of the same or- 
der of magnitude as the anharmonic corrections we have 
calculated for this system. 

Figure [2] shows that convergence with respect to su- 
percell size is reached with smaller cells for H 7 Li and 
D Li than for diamond. In diamond, the phonons with 
the largest anharmonic contribution to the energy are lo- 
cated in the third phonon branch, corresponding to the 
highest energy acoustic modes. These modes represent 
large-scale atomic displacements and larger supercells are 
needed to describe them. In contrast, the largest anhar- 
monic contributions in H 7 Li and D 7 Li arise from the op- 
tical branches, which represent small scale displacements 



of the atoms within the primitive cells. This local an- 
harmonic behavior is reflected by convergence at smaller 
supercell sizes. 



B. Band gap renormalization 

We have used the phonon expectation value formal- 
ism to investigate the ZP renormalization and tempera- 
ture dependence of band gaps in diamond and isotopes 
of lithium hydride. In both cases we have selected the 
thermal band gap (smallest band gap), which is indirect 
for diamond and direct for the lithium hydride systems. 
For diamond, the top of the valence band is at the elec- 
tronic r point, and the bottom of the conduction band 
is along the line joining the T and X points. The ther- 
mal gap is at the electronic X point for the lithium hy- 
dride isotopes. In all calculations reported in this section 
the anharmonicity has been treated at the independent 
phonon level. 

We show the temperature dependence of the electronic 
thermal band gap E g of diamond in Fig. HJ The black 
diamonds are experimental results from ReTT|44l The ZP 
band gap renormalization is E^ F = —462 meV. This 
value is 92 meV larger than that estimated from the 
asymptotic behavior at high temperature in Ref . 05J The 
analysis in Ref. 45 suffers from the lack of experimental 
data beyond the Debye temperature, where the asymp- 
totic behavior should occur, and an analytic functional 
form with the correct asymptotic behavior was fitted to 
the data instead. We attribute some of the discrepancy to 
the inaccuracy of the fitted function used in Ref. 051 Our 
first-principles results can also be compared with other 
theoretical studies of the thermal gap of diamond^ based 
on the empirical pseudopotential methocP^ which find 
a ZP renormalization of approximately —620 meV. The 
calculations reported in Ref. 06 lead to an overestimate of 
the reduction in band gap with temperature compared to 
experiment. This could explain the larger ZP renormal- 
ization found in Ref. 0J5] compared to our first-principles 
results. The band gap renormalization has an additional 
intrinsic contribution from the change in volume of the 
crystal with temperature. For diamond this contribution 
is very small, amounting to a decrease of 16 meV for the 
temperature range in Fig. [4] 

Figure [5] shows the temperature dependence of the 
thermal band gaps of H 7 Li and D 7 Li. The ZP renor- 
malizations arising from electron-phonon interactions are 
Eg F = —84 meV and E g F = —62 meV, respectively; and 
those arising from the change in volume with tempera- 
ture are E g F = +40 meV and E g F = +33 meV, respec- 
tively The isotope effect for the lighter compound leads 
to a larger band gap renormalization for both electron- 
phonon and volume contributions. This is in agreement 
with the expectation that lighter isotopes show larger 
atomic displacements. 

It is interesting to note in the lithium hydride systems 
that the ZP electron-phonon renormalizations are nega- 
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FIG. 4. (Color online) Temperature dependence of the ther- 
mal band gap E s of diamond. The DFT result (red solid 
curve) is offset to match experimental data (black diamonds) 
at zero temperature. The experimental data are from Ref. 1441 



tive (gap closing), while opening of the gap is observed 
at finite temperature. To investigate this we have calcu- 
lated the contribution to the electron-phonon renormal- 
ization from each phonon mode, and the results indicate 
that low-energy phonons tend to open the gap whereas 
high-energy phonons tend to close it. At the ZP level, 
all phonons contribute, and the high energy gap-closing 
phonons dominate. However, as the temperature is in- 
creased, the lower energy phonons are excited more eas- 
ily, leading to opening of the gap. At higher tempera- 
tures, when the high-energy phonons are excited as well, 
the trend is expected to revert to gap closing. In Fig. 
[5| the higher temperatures explored (close to the melt- 
ing point) only show a leveling off of the electron-phonon 
correction. Extending the calculations to higher unphys- 
ical temperatures (not shown) does take the system into 
the gap-closing regime as expected. 

Our first-principles calculations for lithium hydride 
are, as far as we are aware, the first to show non- 
monotonic behavior in the temperature dependence of 
a band gap due to electron-phonon interactions. This 
behavior has been observed experimentally in ternary 
semiconductors^^ anc j it has been described by an an- 
alytic functional form consisting of the sum of two Bose- 
Einstein oscillators of opposite sign and with energies 
related to the low and high frequency phonon modest 
Our calculations show that this difference in behavior be- 
tween the low and high frequency phonons is responsible 
for the nonmonotonic variation of the gap with temper- 
ature. 

In Fig. [5] we also show the sum of the contributions 
to the gap renormalization from electron-phonon inter- 
actions and thermal expansion for H 7 Li and D 7 Li. This 
simple addition of the two terms is the standard proce- 
dure used in the literature.^ However, this approxima- 
tion might not be very accurate for a system such as 
lithium hydride with a large thermal expansion (see Sec. 




400 

Temperature (K) 

FIG. 5. (Color online) Temperature dependence of the ther- 
mal band gap E g of lithium hydride for the isotopes H 7 Li (red) 
and D 7 Li ( green) calculated within DFT. The black double- 
dashed-dotted line indicates the value of the static band gap. 
Dashed lines refer to the renormalized gaps due to electron- 
phonon interactions, dashed-dotted lines refer to the renor- 
malization due to thermal expansion with temperature, and 
the solid lines are the sum of the two contributions. 



VC below). We calculated the electron-phonon inter- 
action renormalization to the band gap at the volume 
including ZP motion for H 7 Li, and the real renormaliza- 
tion is 4 meV larger than the one found by the simple 
addition procedure. 

In Fig. [I| a constant energy shift has been added to the 
theoretical curve to match experiment at zero tempera- 
ture. Standard approximations to the DFT exchange- 
correlation functional such as the LDA or PBE do not 
reproduce the derivative discontinuity with respect to 
particle number of the exact functional. This leads to 
a severe under estim ation of band gaps in semiconductors 
and insulators! 50 * 51 ! This problem can be addressed in var- 
ious ways, the most straightforward being the application 
of a scissor operatoj 50 * 52 ! to the band gaps. When calcu- 
lating band-gap renormalization, we are interested in the 
change in the band gap, rather than its absolute value. 
The fact that all atomic configurations suffer from a 
similar gap underestimation suggests that the calculated 
changes may be accurate even when using standard func- 
tionals. This is supported by the calculations of Giustino 
and co-workers^ of the electron-phonon band-gap renor- 
malization in di amon d using Allen-Heine-Cardona per- 
turbation theoryP^EH They use both the LDA and the 
LDA corrected with a scissor operator, obtaining consis- 
tent results. Note that Giustino and co-workers look at 
the direct optical gap instead of the therm al gap that we 
have studied. Cannuccia and MarinpMSiJ U se many-body 
perturbation theory to study the optical gap of diamond 
as well, and they report a significant structure in the 
spectral function beyond the quasiparticle picture, corre- 
sponding to subgap polaronic states. 
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FIG. 6. (Color online) Temperature dependence of the coeffi- 
cient of linear expansion a for diamond. The black diamonds 
are experimental results from Ref. 1551 The inset shows the 
temperature dependence of the lattice parameter. 



FIG. 7. (Color online) Temperature dependence of the coef- 
ficient of linear expansion a for lithium hydride with isotopic 
compositions H 7 Li and D 7 Li. The experimental data are from 
Ref. 1561 The inset shows the temperature dependence of the 
lattice parameters. 



C. Thermal expansion 



We have used the formalism described in Sec. IIII CI to 
calculate the internal vibrational stress and the corre- 
sponding equilibrium configuration of the crystal as a 
function of temperature. The structures of diamond, 
H 7 Li, and D 7 Li are described by a single parameter, 
so the different equilibrium configurations differ only in 
the volume. We present results for the thermal expan- 
sion of diamond, H 7 Li, and D 7 Li. All calculations re- 
ported in this section have treated anharmonicity at the 
independent-phonon level. 

In Fig. [6] we show the temperature dependence of the 
lattice parameter a{T) of diamond and the correspond- 
ing temperature dependence of the linear coefficient of 
thermal expansion, 



a(T) = 



1 da 
adT' 



(25) 



We compare our results with experimental data from Ref. 
|5"51 The agreement is good for temperatures below 1200 
K, but at higher temperatures the calculated a{T) under- 
estimates the experimental values. This effect is also seen 
in calculations using the quasiharmonic approximation.^ 
Our methodology assumes a volume-independent inter- 
nal vibrational stress tensor for the volumes of interest 



at each iteration, as described in Sec. HI C Iterating the 
calculation at the volume including the ZP phonon pres- 
sure, leads to a change in the phonon pressure of dia- 
mond of less than —0.05 GPa, which represents less than 
1.5% of the total phonon pressure. This means that for 
diamond a single iteration is sufficient for obtaining con- 
verged results. 

The results for the lithium hydride isotopes are shown 
in Fig. [7] and are compared with experimental data from 
Ref. [5J3 The lattice parameters of the different isotopes 



show the expected behavior, with the lighter compound 
leading to a larger ZP phonon pressure and therefore to 
a larger lattice expansion. In contrast, the linear coeffi- 
cient of thermal expansion is larger for the heavier deu- 
terium compound, in agreement with experiment. The 
phonons of D 7 Li have lower energies than those of H 7 Li, 
and therefore they are excited more easily with increas- 
ing temperature. This leads to the thermal expansion 
coefficient being larger for D 7 Li at low temperatures. 

The volume change in lithium hydride due to ZP mo- 
tion is much larger than that in diamond. The change 
in ZP phonon pressure between the first and second it- 
erations is as large as —0.22 GPa, representing 11% of 
the absolute value. In this case the second iteration is 
important for obtaining accurate results. The change 
in pressure is negative because the expanded lattice (af- 
ter the first iteration) has softer phonons than the static 
DFT lattice, which leads to smaller phonon pressures. 

The absolute values of the lattice parameters within 
DFT are known to deviate systematically by a few per- 
cent from experimental results. However, the differences 
in lattice parameters within a given DFT approximation 
are expected to be accurate. These are displayed in Ta- 
ble |Hj Mounet and Marzari (Ref. find the same ZP 
expansion for diamond using the PBE functional, rather 
than the LDA that we have used, even though the abso- 



TABLE II. Static DFT and ZP renormalized lattice parame- 
ters of diamond, H 7 Li, and D 7 Li. 





a static ( a u _) 


a zp (a.u 


) Aa (a.u.) 


Diamond (LDA) 


6.669 


6.694 


+0.025 


H 7 Li (PBE) 


7.600 


7.758 


+0.158 


D 7 Li (PBE) 


7.600 


7.726 


+0.126 
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lute value of the lattice parameter is somewhat different. 
Anderson and co-workers^ report experimental results 
for the lattice parameter of H 7 Li and D 7 Li at a temper- 
ature of 83 K, finding a difference between the value for 
the two isotopes of 0.032 a.u., in good agreement with 
the value of 0.034 a.u. we have calculated. 



VI. CONCLUSIONS 

We have described a formalism for studying tempera- 
ture dependent anharmonic vibrational properties in pe- 
riodic systems using first-principles calculations. The BO 
energy surface is mapped using a principal axes approx- 
imation, and the resulting equations are solved using a 
Hartree mean-field approach and second-order perturba- 
tion theory. We then use the vibrational anharmonic 
wave function to calculate expectation values of phonon- 
dependent quantities. In particular, we have studied the 
effects of electron-phonon interactions on electronic band 
gaps, and the role of the stress tensor in thermal expan- 
sion. Other quantities can be studied within this frame- 
work, such as the mean atomic positions within the unit 
cell. 

The methodology we have presented involves the solu- 
tion of the electronic Schrodinger equation at a set of low 
energy atomic configurations. Our results are based on 
density functional theory calculations, but the methodol- 
ogy is independent of the particular electronic structure 
method used. For instance, for larger systems it might 



be necessary to use force fields, or, if higher accuracy is 
required, quantum Monte Carlo methods!^ 

We have tested the formalism on diamond, lithium hy- 
dride and lithium deuteride, which exhibit quite small 
anharmonicities. However, while diamond is nearly har- 
monic, the small anharmonicity in lithium hydride arises 
from a cancellation of the contributions from single- 
phonon and two-phonon terms. The results for both 
band-gap renormalization and lattice expansion as a 
function of temperature have been found to be in good 
agreement with experimental results for diamond. For 
H 7 Li and D 7 Li, our thermal expansion results also agree 
with experiment but, as far as we are aware, there are 
no experimental data for the temperature dependence of 
the band gap and our results serve as a prediction that 
could be tested. 

In our calculations using the VSCF method plus per- 
turbation theory we appear to have converged the results 
with respect to the description of anharmonicity. This 
suggests that our approach can be used to estimate an- 
harmonic effects in systems with substantially stronger 
anharmonicity than diamond and lithium hydride. 
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